#ifndef LMMPRICER_HPP
#define LMMPRICER_HPP

#include "swap.hpp"
#include "BermudanSwaption.hpp"
#include "LMMAdvancer.hpp"
#include "LMMPathPool.hpp"
#include "leastSquareSolver.hpp"
#include "BlackScholesPricer.hpp"
#include "varianceStructure.hpp"
namespace TermStructure
{
	template<typename TVarianceStructure,template<typename, typename,typename> class TMultiDimensionGaussian, typename T = double,typename TSize = long>
	class LMMPricer
	{
	protected:
		TimeSequence<T,TSize> *_ts;
		VarianceStructure<T> *_vs;
		LMMPathPool<TVarianceStructure, TMultiDimensionGaussian, T, TSize> _pathpool;
	public:
		LMMAdvancer<TVarianceStructure,TMultiDimensionGaussian, T,TSize> UnderlyingAdvancer;
		void SetLMMAdvancer( TimeSequence<T,TSize> &ts, TVarianceStructure &variance)
		{
			this->UnderlyingAdvancer.InitiateForwardRateSequence();
			this->UnderlyingAdvancer.Set(variance, ts);
			this->_ts=&ts;
			this->_vs=&variance;
		}
		ForwardRateSequence<T,TSize> &InitialForwardRateSequence()
		{
			return UnderlyingAdvancer.CurrentSequence();
		}
		const ForwardRateSequence<T,TSize> &InitialForwardRateSequence() const
		{
			return UnderlyingAdvancer.CurrentSequence();
		}
		T PriceEuropeanSwaption(const Swap<T,TSize> &swap, const T &strike, TSize chooserIndex, TSize N = 1000000)
		{
			T result = 0;
			UnderlyingAdvancer.AnchorIndex = chooserIndex;
			UnderlyingAdvancer.CurrentTime = 0;
			UnderlyingAdvancer.TimeStep = UnderlyingAdvancer.GetTimeSequence().GetTime(chooserIndex);
			UnderlyingAdvancer.PreAdvance();
			swap.Strike = strike;
			for(TSize i = 0; i < N; i++)
			{
				UnderlyingAdvancer.SingleAdvance();
				T v = swap.Price(UnderlyingAdvancer.NextSequence());
				if(v >= 0) result += v/N;
			}
			UnderlyingAdvancer.CurrentSequence().EstablishDiscountFactors();
			return result * UnderlyingAdvancer.CurrentSequence().DiscountFactorAt(chooserIndex);
		}

		T PriceTriggerSwap(const Swap<T,TSize> &swap, const T &strike,TimeSequence<T,TSize> &chooserTime,TSize N = 1000000)
		{
			T result = 0;
			T v =0;
//			UnderlyingAdvancer.AnchorIndex = chooserIndex;
			swap.Strike = strike;
			
			ForwardRateSequence<T,TSize> **fwd_path = new ForwardRateSequence<T,TSize>*[N] ;
			for (TSize n = 0; n < N;n++){
				fwd_path[n] =  new ForwardRateSequence<T,TSize> [chooserTime.GetSize()+1];
			}
			
			UnderlyingAdvancer.CurrentTime = 0;
			
			for (int i = 1; i < chooserTime.GetSize(); i++)
			{
				UnderlyingAdvancer.TimeStep = chooserTime.GetTime(i);;
				UnderlyingAdvancer.PreAdvance();
				for (TSize n = 0; n < N;n++){
					UnderlyingAdvancer.SingleAdvance();
					fwd_path[n][i] = UnderlyingAdvancer.NextSequence();
				}
			}

			for(TSize i = 0; i < N; i++)
			{
				for (int j = 1; j < chooserTime.GetSize(); j++)
				{
					
					T v = swap.Price(fwd_path[i][j]);
					
					if(v >= 0){
						UnderlyingAdvancer.CurrentSequence().EstablishDiscountFactors();
						result += v/N/(UnderlyingAdvancer.CurrentSequence().DiscountFactorAt(UnderlyingAdvancer.AnchorIndex,UnderlyingAdvancer.CurrentTime));
						break;
					}
				}
			}
			UnderlyingAdvancer.CurrentSequence().EstablishDiscountFactors();
			return result * UnderlyingAdvancer.CurrentSequence().DiscountFactorAt(UnderlyingAdvancer.AnchorIndex);
		}

		T PriceConstMaturitySwap(const T &fwp, Swap<T,TSize> &swap, const TSize swap_b, const T &strike ,TSize N = 50000)
		{
			T s =0;
			T v =0;  
			swap.SwapEndIndex = swap_b;
			swap.SwapEnterIndex = fwp;
			swap.Strike =strike ;
				
			UnderlyingAdvancer.CurrentTime = 0;
			ForwardRateSequence<T,TSize> **fwd_path = new ForwardRateSequence<T,TSize>*[N] ;
			for (TSize n = 0; n < N;n++){
				fwd_path[n] =  new ForwardRateSequence<T,TSize> [fwp+1];
			}

			for (TSize i =1; i <= fwp;i++){
				 
				UnderlyingAdvancer.TimeStep = _ts->GetTime(i);
				UnderlyingAdvancer.PreAdvance();
				for (TSize n = 0; n < N;n++){
					UnderlyingAdvancer.SingleAdvance();
					fwd_path[n][i] = UnderlyingAdvancer.NextSequence();
				}
			}

			for (TSize n = 0; n < N;n++){

			for (TSize i =1; i <= fwp;i++)
			{
				swap.SwapEndIndex = swap_b+i -1;
				swap.SwapEnterIndex = fwp+i -1;
				s = swap.Price(fwd_path[n][i]);
				v = fwd_path[n][i][i];
					
				for (TSize j = i+1;j<=fwp;j++)
				{
					s = s*(1+fwd_path[n][int(fwp)][j]*UnderlyingAdvancer.GetTimeSequence().TimeDifference(j));
				}			
				v += (v - s)/N;
			}
			
			}
			UnderlyingAdvancer.CurrentSequence().EstablishDiscountFactors();
			return v*UnderlyingAdvancer.CurrentSequence().DiscountFactorAt(fwp);

		}

		T PriceRatchetcaps(const T &fwp, TSize N = 50000)
		{
			T v =0;
			ForwardRateSequence<T,TSize> *fwd_path_end = new ForwardRateSequence<T,TSize>[N];
			UnderlyingAdvancer.CurrentTime = 0;
			UnderlyingAdvancer.TimeStep = UnderlyingAdvancer.GetTimeSequence().GetTime(fwp);
			UnderlyingAdvancer.PreAdvance();

			for (TSize n = 0; n < N;n++){
				UnderlyingAdvancer.SingleAdvance();
				fwd_path_end[n] = UnderlyingAdvancer.NextSequence();
			}

			for (TSize n = 0; n < N;n++){
				UnderlyingAdvancer.TimeStep = _ts->TimeDifference(1);
				UnderlyingAdvancer.Advance();
			for (TSize i =0; i <= fwp-3;i++)
			{
			   
			 	v = (UnderlyingAdvancer.NextSequence()[i+2]-UnderlyingAdvancer.CurrentSequence()[i+1]);
				v = max<T>(v,0);
				if (v >0)
				{
					for (TSize j = i+2;j<=fwp;j++)
					{
						v = v*(1+fwd_path_end[n][j]);				}			
					}
				v += (v)/N;
				}
			}
			UnderlyingAdvancer.CurrentSequence().EstablishDiscountFactors();
			return v*UnderlyingAdvancer.CurrentSequence().DiscountFactorAt(fwp);

		}
		
		T PriceCaplet(T time_now,T time_beg,T time_end,T K,TSize N = 50000)
		{
			T v = 0;
			TSize end_i = _ts->GetIndex(time_end);
			TSize beg_i = _ts->GetIndex(time_beg);
			UnderlyingAdvancer.CurrentTime = time_now;
			UnderlyingAdvancer.TimeStep = time_beg;
			UnderlyingAdvancer.PreAdvance();
			UnderlyingAdvancer.SingleAdvance();
				
			for(int i = 0; i < N; i++)
		    {
		    UnderlyingAdvancer.SingleAdvance();
			 if(UnderlyingAdvancer.NextSequence()[end_i] >= K) {
				UnderlyingAdvancer.NextSequence().EstablishDiscountFactors();
				v += _ts->TimeDifference(end_i)*(UnderlyingAdvancer.NextSequence()[end_i] - K)/(1+UnderlyingAdvancer.NextSequence()[end_i]*_ts->TimeDifference(end_i))/N/(UnderlyingAdvancer.NextSequence().DiscountFactorAt(end_i,beg_i));
			 }
		 }
			UnderlyingAdvancer.CurrentSequence().EstablishDiscountFactors();
			return v * UnderlyingAdvancer.CurrentSequence().DiscountFactorAt(end_i);
		}
		T PriceSwap(Swap<T,TSize> &swap,const T time_now, const TSize swap_b,const TSize swap_a,const T strike,TSize N = 50000)
		{
			T s = 0;
			swap.SwapEndIndex = swap_b;
			swap.SwapEnterIndex = swap_a;
			swap.Strike =strike ;
			UnderlyingAdvancer.CurrentTime = time_now;
			UnderlyingAdvancer.TimeStep = _ts->GetTime(swap_a)-time_now;
			UnderlyingAdvancer.PreAdvance();
			UnderlyingAdvancer.SingleAdvance();
			s = swap.Price(UnderlyingAdvancer.NextSequence());

			return s;
		}



		public:
		T PriceBermudanSwaption(Swap<T,TSize> &contract, TSize NumberOfPath = 10000)
		{
			T result = 0;
			this->_pathpool.SetLMMAdvancer(this->UnderlyingAdvancer);
			this->UnderlyingAdvancer.AnchorIndex = contract.SwapEndIndex;
			this->_pathpool.NumberOfPath = NumberOfPath;
			this->_pathpool.CurrentTimeIndex = 0;
			this->_pathpool.EndIndex = contract.SwapEndIndex - 1;
			this->_pathpool.Generate();
			SimpleMatrix<T,TSize> cvx(NumberOfPath, 6);
			SimpleMatrix<T,TSize> cvy(NumberOfPath,2);
			SimpleMatrix<T,TSize> coef(6,1);
			Swap<T,TSize> anotherContract(this->UnderlyingAdvancer.GetTimeSequence(), contract.SwapEndIndex, 0);
			anotherContract.Strike = contract.Strike;
			anotherContract.SwapEnterIndex = contract.SwapEndIndex  - 1;
			for(TSize i = 0; i < NumberOfPath; i++)
			{
				cvy(i,0) = anotherContract.Price(this->_pathpool(i, contract.SwapEndIndex - 1), contract.SwapEndIndex - 1);
				if(cvy(i,0)<0) cvy(i,0) = 0;
				else cvy(i,0) /= this->_pathpool(i,contract.SwapEndIndex - 1).DiscountFactorAt(contract.SwapEndIndex, contract.SwapEndIndex - 1);

			}
			
			for(TSize i = contract.SwapEndIndex - 2; i > contract.SwapEnterIndex; i--)
			{
				//std::ofstream out("tester.csv");
				anotherContract.SwapEnterIndex = i;
				for(TSize pi = 0; pi < NumberOfPath; pi++)
				{
					cvy(pi,1) = anotherContract.Price(this->_pathpool(pi,i),i);
					cvy(pi,0) *= this->_pathpool(pi,i).DiscountFactorAt(contract.SwapEndIndex, i);
					T x = _pathpool(pi,i).SwapRate(i+1, contract.SwapEndIndex);
					T y = _pathpool(pi,i)[i+1];
					cvx(pi,0) = 1;
					cvx(pi,1) = x;
					cvx(pi,2) = y;
					cvx(pi,3) = x*y;
					cvx(pi,4) = x*x;
					cvx(pi,5) = y*y;
					//out<<x<<","<<y<<","<<cvy(pi,0)<<std::endl;
				}
				LeastSquareSolver<T,TSize>::Fit(cvx,cvy,coef);
				//outputTo(cout,coef);
				for(TSize pi = 0; pi < NumberOfPath; pi++)
				{
					T temp = 0;
					for(TSize l = 0; l < 6; l++)
					{
						temp += cvx(pi,l)*coef(l,0);
					}
					if(temp < 0) temp = 0;
					cvy(pi,0) = temp > cvy(pi,1) ? temp: cvy(pi,1);
					cvy(pi,0)/= this->_pathpool(pi,i).DiscountFactorAt(contract.SwapEndIndex, i);
				}
				//out.close();
			}
			for(TSize pi = 0; pi < NumberOfPath; pi++)
			{
				result += cvy(pi,0)/NumberOfPath;
			}
			_pathpool(0,0).EstablishDiscountFactors();
			return result*_pathpool(0,0).DiscountFactorAt(contract.SwapEndIndex, 0);
		}

		T PriceCaption(const T &cap_end, const T &cap_strike, const T &strike, TSize chooserIndex, TSize N = 1000000)
		{
			T v = 0;
			UnderlyingAdvancer.CurrentTime = 0;
			UnderlyingAdvancer.TimeStep = UnderlyingAdvancer.GetTimeSequence().GetTime(chooserIndex);
			UnderlyingAdvancer.PreAdvance();
			UnderlyingAdvancer.AnchorIndex = chooserIndex;
			
		    BlackScholesPricer<VarianceStructure<T,TSize>, T> bsp;
			bsp.SetTimeSequence(*_ts);
			UnderlyingAdvancer.CurrentSequence().EstablishDiscountFactors();
			for (int n = 0; n < N ; n ++) {
				UnderlyingAdvancer.SingleAdvance();
				v = max<T>(bsp.PriceCap(UnderlyingAdvancer.NextSequence(), *_vs, cap_strike,cap_end,chooserIndex, chooserIndex)-strike,0);
				v += v/N;	
			}

			return v* UnderlyingAdvancer.CurrentSequence().DiscountFactorAt(UnderlyingAdvancer.AnchorIndex);
		
		
		}
		


	};
} 

#endif


/*
	TimeSequence<double> ts;
	for(int i = 1; i <= 40; i++)
	{
		ts.PutTime(i*0.25);
	}

	ForwardRateSequence<double> frs;
	frs.SetTimeSequence(ts);
	for(int i = 0; i < ts.GetSize(); i++)
	{
		frs[i] = 0.05;
	}

	VarianceStructure<double> vs;
	vs.d=0.05;
	vs.a=0.09;
	vs.b=0.44;
	vs.c=0.11;
	vs.beta=0.1;
	Swap<double> swap(ts, 40, 4);


	LMMPricer<VarianceStructure<double>,PCABasedMultiDimensionGaussian,double> pricer; //PCABasedMultiDimensionGaussian
	pricer.SetLMMAdvancer(ts, vs);

	for(int i = 0; i < pricer.InitialForwardRateSequence().GetSize(); i++)
	{
		pricer.InitialForwardRateSequence()[i] = 0.05;
	}
//	cout<<"European Swaption:"<<pricer.PriceEuropeanSwaption(swap, 0.05, 4)<<endl;
	TimeSequence<double> chooser;
	for (int i =1;i<4;i++)
	{
		chooser.PutTime(i*0.25);
	}
	Swap<double> swap1(ts);

//	cout<<"Trigger Swaption:"<<pricer.PriceTriggerSwap(swap, 0.05, chooser,5000)<<endl;
//	cout<<"Const Maturity Swap"<<pricer.PriceConstMaturitySwap(10,swap1, 22,2,0,1000)<<endl;
//	cout<<"Ratchet caps"<<pricer.PriceRatchetcaps(10,1000)<<endl;
//	cout<<"caplet:"<<pricer.PriceCaplet(0,0.25,0.5,0.05,100000)<<endl;
	cout<<"swap:"<<pricer.PriceSwap(swap,0,40,4,0.05,100000)<<endl; //swap(0;1,2.5;K=0)
*/